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ABSTRACT 



Continuous and discrete wavelet transforms are discussed. Recursive generations of 
compactly supported wavelet time functions are described. Haar, Daubechies, and complex 
modulated Gaussian wavelets are utilized to analyze phase shifts in signals with noise. 
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1 . INTRODUCTION 

We have investigated several time-frequency decomposition techniques and their 
applications in the detection and analysis of digitally modulated signals [10], [14], [13]. In this 
report, we will specifically discuss wavelet decompositions and present a few preliminary results. 

The one particular digital modulation that we consider in here is binary phase-shift-key 
(BPSK) signal, which is a sinusoid modulated by a baseband rectangular data waveform. Over a 
duration of one binary data bit, the phase of a BPSK wave is either 0' or 180“ depending on the 
value of the data, which is either +1 or -1 . For example, a wave with data rate 9600 bps can 
change phase every 1/9600 seconds. In tactical communications using direct-sequence spread- 
spectrum (DS/SS) [1 1], the corresponding spread BPSK wave can change phase N times faster, 
where N is the spectrum spreading factor. Since BPSK and its extensions (such as QPSK and 
OQPSK) are the most predominant modulation formats used in practice, there are interests in 
developing efficient methods for detecting the phase shifts under various channel and noise 
scenarios. Previous research efforts have used cyclostationarity [6] and ambiguity function with 
specific kernel as the smoothing function [23], [13]. 

In Chapter 2, we will discuss the so-called continuous wavelet transform (CWT). In 
Chapter 3, we will discuss the discrete wavelet transform (DWT) and the computations of 
orthonormal compactly supported scaling and wavelet functions on the dyadic numbers. In 
Chapter 4, we use the wavelet time-frequency decomposition techniques to detect ±180“ phase 
shifts in noise. A few concluding remarks are given in Chapter 5. 
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2 . CONTINUOUS WAVELET TRANSFORM 

The continuous wavelet transform (CWT) of a L2(E) signal x(t) is defined as [8], [2], [3]: 

X(t, a) = -^ [ x(u) 1 du 

Va J-oo \ a / (2.1) 



where \}i(t) is the wavelet function, " * ” denotes complex conjugation, and a > 0. The scaling 
factor "a" in (2.1) is taken to be power of 2 in this report, i.e., a = 2~^, for m = 0, 1, 2, 3 , . . . 
The inverse transform is given by: 





(2.2) 


where the constant 


(2.3) 


must be finite. This implies that the wavelet function must satisfy 




^(0 = 0, i.e., J \K(t)dt =0 


(2.4) 



Thus, the wavelet function is bandpass type with Fourier transform 4^(f) = 0 at f = 0 ; 
equivalently, the area under the wavelet function is exactly 0. Note that we can rewrite (2.1) as 

X(t,a) = VaJ X(f) 'F*(af) exp(j27ift) df (2.5) j 

in terms of the Fourier transforms. The CWT can be computed via (2.1) or (2.5). j 

t 

The CWT X(t, a) represents a decomposition of x(t) into a time-scale plane. Usually, 

X(t,a) is graphed against time t and the log scale log 2 (a). The CWT can be regarded as a constant- 
Q bandpass filter since (2.1) is a convolution of two time functions. The linear bandpass filter i ’ 

impulse response at level m is ' 

Vm(t) = 2"^ V*(2"'t) (2.6) ' I 

If we define "fm" as the r.m.s. center frequency of the passband by 

oo 

f (f-fJ^I'Fjf)Pdf =0 

Jo 



(2.7) 
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then we have 



fm = 2^ fo and Om = 2<^ Oq 



( 2 . 8 ) 



where fo is the r.m.s. center frequency and 2oo is the r.m.s. bandpass bandwidth of the basic 
wavelet function \|i(t). The ratio of fn, to Om is always equal to 



which is constant for all decomposition level m. As the level of decomposition (i.e., m) increases, 
the filtering function ym(0 of (2.6) becomes narrower in time. The time-resolution of the wavelet 
decomposition therefore increases with m. On the other hand, the r.m.s. bandpass bandwidth 
(i.e., 20m ) of \ltm(t) increases exponentially with m, and therefore the frequency-resolution of the 
wavelet decomposition is better for smaller m. Unlike the classical time-frequency transforms 
such as the short-time Fourier transform (STFT) or the Gabor transform (GT) [1], [9], which do 
not scale the filtering function and and have constant time-resolution at all frequencies, the wavelet 
transform has higher time-resolution at higher frequency. One can observe that for CWT, its 
time-resolution is increasingly better at higher frequencies, while its frequency-resolution is 
increasingly better at lower frequencies. 

We will use the complex modulated Gaussian wavelet function; 



Of course the Fourier transform of a Gaussian function is a Gaussian function, and the complex 
exponential time function exp(j27tfyt) will cause a shift by in frequency in the Fourier 
transform. It is then possible to choose fy big enough so that for complex modulated Gaussian 
wavelet we have 



example, let x(t) = cos 27ifct, for t in a very large, yet finite, window. Using (2.5), we have 




(2.9) 



\|/(t) = exp(j27ifv|it) exp(-t2/2o2) 

= cos 27tfvj>t exp(-t2/2a^) + j sin 27tfyt t\p(-tV2o^) 



( 2 . 10 ) 



'F(f) = 0 for f < 0 



( 2 . 11 ) 



Condition (2.1 1) is especially desirable when the wavelet is used to analyze sinusoids. For 



X(t, a) = ^ exp(j27tfc0 T*(afc) + ^ exp(-j27tfct) 'F*(-afc) 



( 2 . 12 ) 
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If 'F(f) satisfy (2.1 1), then the second term in (2.12), which is a contribution by the impulse 
6(f+fc), drops out. It follows that the phase of the resulting CWT X(t, a) will nicely reveal the 
center frequency of the carrier wave x(t) = cos 2jrfct. For other stationary or non-stationary 
signals, magnitude and phase of X(t, a), with V|/(t) being the complex Gaussian wavelet, will reveal 
the frequency as well as the time behaviors of the signal. 

We also remark that we can multiply a suitable baseband function b(t) with the complex 
sinusoid exp(j2jtfyt) to obtain a wavelet function: 



so that its Fourier transform is zero on the negative frequency axis. The above construction yields 
a large class of potentially useful wavelet functions. The preference to use the Gaussian envelope 
exp(-t2/2o2) is due to the fact that the Gaussian function achieves the lower bound in the 
uncertainty principle: 



with equality, where At is the r.m.s. time-spread and Af is the r.m.s. frequency-spread of the 
signal [20]. Note that two pulses in time can be discriminated only if they are more than At 
seconds apart in time, and two sinusoids can be discriminated only if they are more than Af apart 
in frequency. Equation (2. 14) states that one can only trade time-resolution for frequency- 
resolution, or vice versa. 

Using CWT, the frequency-resolution Af is in the order of Om = 2^ Gq. Viewing the 
CWT X(t, a) in a time-scale plane allow us to look at the signal at different scale (i.e., different 
frequencies in the passband of 2am Hz wide around the center frequency fm). The CWT indeed 
presents a multiresolutional analysis of the signal, as if the signal is being processed by a bank of 
filters with logarithmically spaced centered frequency fm but constant bandwidth to frequency ratio 
of 2am/ffn = 2ao/fo. These relationships are best visualized by considering the complex modulated 
sine wavelet function: 



\|t(t) = exp(j27tfv|,t) b(t) 



(2.13) 



AtAf > — 
4k 



(2.14) 




sinc^ 



(2.15) 
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which has a rectangular passband with bandwidth equal to 2cq = W. The center frequency of V|/(t) 
is exactly equal to fy. We also propose to use this complex modulated sine wavelet to analyze our 
signals. 

Computation of X(t, a) could be lengthy and expensive. In practice, however, we only 

have a finite-length sampled sequence: x(nT), n = 0, 1, 2, Ng-l, where Ns is the 

sequence length and the samphng rate is at least equal to the Nyquist rate, which is twice the 
highest frequency component in the signal x(t). We will only compute X(t, a) at the sampling 
times t = nT and at levels m = 0, 1, 2, . . . . , [log 2 Ns]. The integral of (2.1) is now reduced to a 
finite summation and the computation burden is greatly relieved. In this discrete-time 
approximation of the CWT, we essentially perform a multiresolutional analysis of an auxiliary 
signal Xa(t), which is a piecewise linear function constructed by joining two consecutive sampled 
points x(nT) and x((n+l)T) with a straight line. Towards this end, the CWT of (2.1) produces a 
total of Nmax = [log2 Ns] + 1 snap-shots of the auxiliary signal via bandpass filterings. 

The complex modulated Gaussian wavelet of (2.10) is a nice analytical function. It is 
continuous and differentiable everywhere. Perhaps one drawback is that it has an infinitely span 
and therefore it is not a causal filter impulse response. In practice, however, this mishap can be 
resolved by using a wavelet with sufficiently fast decay over the processing window so that it can 
be regarded as of finite span. We select the Gaussian wavelet because of two additional reasons: 
(a) tunability; and (b) bandwidth adjustability. By varying the modulating frequency fy, we can 
place the passband anywhere on the frequency axis. Second, by altering o, the bandwidth of the 
wavelet function \|r(t) can be selected at will. Each (fy, a) selection will result in a 
multiresolutional analysis that consists of Nmax snap-shots of the signal. A systematic selection 
of (fy, a) will allow a thorough analysis of the signal at all possible time-resolution and 
frequency-resolution levels. 

We show an example of a complex modulated Gaussian wavelet function \j/(t) with fy = 8 
Hz and o = 0.12 in Figure 2.1. This value is selected so that all the energy in the wavelet is 
essentially contained in the compact time interval [-0.5, 0.5] for presentation purpose. In the left 
column of Figure 2.1, the real part, imaginary part, magnitude, and phase of the wavelet are 
depicted. In the right column, the energy spectral densities (which is magnitude square of the 
Fourier transform) of the real part, imaginary part, and the wavelet function itself are shown. Note 
that this complex wavelet has no spectral component in the negative frequency axis, a desirable 
property as discussed previously. The spectrum is centered around f = fy = 8 and the r.m.s. 
passband bandwidth of \|/(t) is 2a = 0.24. In processing a given signal x(t), the signal is first 



nonnalized to be contained in the time interval [0, 1]; it is then decomposed by the analyzing 
functions V|/m(t), with t}/o(t) = \jr(t) contained also in the interval [0, 1]. 
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Fig. 2.1 



Complex modulated Gaussian wavelet function \y(t): time-domain functions and 
energy spectral densities. f\(, = 8, a = 0.12. 



3. DISCRETE WAVELET TRANSFORM 



In this section we will discuss an orthonormal decomposition of a signal x(t) using a 
scaling/wavelet function pair, and describe the so-called discrete wavelet transform (DWT) 
algorithm using orthonormal wavelet functions [4], [16], [17], [18]. We will also discuss the 
multiple-phase DWT, and generations of wavelet functions that have compact (i.e., closed and 
bounded) support. 

3.1 Orthonormal Wavelet Decomposition 

Consider a nested sequence of closed subspaces { Vm’, me 2 } in L2(IR) such that 

{0} = - -c V_ 2 C V_, c Voc V, c Vjc- • =L^(1R) (3.1) 

where O = {0} and 0 = L (IR). One can find an unique scaling 

m€ Z m€ Z 

(or averaging) function (J)(t) such that {<|)n,j,(t); n€ 2} 
is an orthonormal basis of Vm where 



= a'"t-nT) 



(3.2) 



by normalizing, dilating or contracting, and time-shifting (j)(t). The value of T (>0) is arbitrary. As 
is common in dyadic multiresolution analysis, we take a (9^0,1) to be 2 in this report. 

The orthogonal projection of a finite-energy deterministic function, x(t), onto the subspace 
Vm is an approximation of that function at resolution level m. We can write: 



~ ^ni,n 
neZ 

where 

t>ni.n = < X(t), <|)m.n(t) > = f Cn W = ( X ♦ q )(n2'"Y) 

(3.3) 



is a coefficient equal to the orthogonal projection of x(t) onto the basis function (]>ni,n(t)- Here, x(t) 
= x(-t). Since [Vm) is a sequence of densed and nested subspaces, we have Xm(t) x(t) in 



1 1 



L2(E) 35 m — > oo. By dilating (when m<0) or contracting (m>0) the scaling function, 
approximation for x(t) at any resolution level can be obtained by increasing the number of 
coefficients exponentially with m. Although the set of functions {<|)m,n(0; (in, n) e spans 
the entire space L2(]R), it is not a minimal set because the subspaces { Vm} are nested; in 
particular, <|)m,n(0 and <j>i,ic(t) are correlated for Jt. For many applications, for example in 
image processing and compact coding, it is desirable to construct an orthonormal multiresolution 
basis for L2(]R). Towards this end, let Vm = Vm-i ® Wm-i, where Wm-i is the orthogonal 
complement of Vm-i in Vm- Consider Vq with {<t>(t-nT); n e Z } as its orthonormal basis. One 
can construct a wavelet function \|/(t) from <j>(t) so that {\jr(t-nT); n e Z } is an orthonormal basis 
for Wq; furthermore, {\|im,n(t); n e Z } = { 2"'^\|/(2“'t - nT); n e Z } is an orthonormal basis for 
Wm, and the entire set { Vm.nCO*. (ni,n) e Z^} is an orthonormal multiresolution basis for L2(E). 
The relationship between the subspaces Vs and Ws is depicted in Fig. 3.1. 




{ 0 } 






Wm-l 



L^(R) 



Figure 3.1 Subspaces Vs and Ws. 



The orthogonal projection of x(t) onto the subspace Wm is the detail function at resolution 
level m; 

^m(^) “ ^m,n M^m,n(^^» 



= X 



m+1 



(t) -xjt). 



~ ^ x(0, Vm,n(^^ ^ 



(3.4) 



This detail function constitutes the difference between the approximation function of x(t) at 
resolution level m+1 and resolution level m. Thus, the coefficients {bm,n’. n e 2 } capture the 
information about x(t) at resolution level m; while the coefficients {dm,n; n e Z } capture the detail 

information going from resolution level m to resolution level m+1. Note that L2(R) = 

and the Ws are mutually orthogonal subspaces. The sequences {bm,n; n e 2 } and {dm,n; n e 2 } 

can be obtained directly through linear filtering of x(t) and then sample at time t = n2‘™T (Fig. 

* * 

3.2). Note that JF [(f)^ ^C-t)] yields a lowpass transfer function, while bandpass. 

♦ 

The wavelet decomposition filters ^(-t) can be interpreted as a parallel bank of bandpass filters 



(with overlapping passbands for adjacent resolution levels) that have constant size passband on the 
logarithm scale. 

t = n2“*^T 







Figure 3.2 Orthonormal wavelet decomposition of the signal x(t). 



3.2 Discrete Wavelet Transform 



There exist recursive (pyramidal) algorithms that efficiently compute the lower resolution 
level coefficients {bi,n‘. n e Z } and {d^,n’» n e 2 } from the higher resolution level coefficients 
{bm.nl n e Z }, i < m; conversely, the bs at the higher resolution level m can be computed from 
the bs at a base (lowest) resolution level M and all the ds at lower resolution levels i, M < i < m 
These pyramidal algorithms are the so-called discrete wavelet transform (DWT) and they are 
described by the block diagrams in Fig. 3.3. 



(a) resolution 



m,n 




2x-downsample decimators 
(keep every other samples) 



(b) reconstruction 




2x-upsample interpolators 
(insert "zero" between each sample) 



Figure 3.3 Discrete wavelet transform. 



In Fig. 3.3, h(k) is the impulse response of a discrete-time conjugate filter that satisfies 
IH((oT) |2 + I H(toT-t-7i) |2 = 1 



(3.5) 



1 



and g(k) is the impulse response of the mirror filter of H with 
G((d) = exp(-jm) H*(o>+7t) 



(3.6a) 



and H(cdT)G*(©T) + H(©T + 7t)G*(coT + n) = 0 



(3.6b) 



One common choice of g(n) that satisfies (3.6) is g(n) = (-I)’"" h(l-n). The scaling function and 
the wavelet function are orthogonal to each other and they are related by 



<J)(2co) = H(co) <I>(co) 
4^(203) = G((0) 0(03) 



(3.7) 



It follows that 0(03) = n^j H(2“"^03); depending on the choice of H(03), it is possible to obtain 



scaling function <|)(t) that has good localization properties in both the frequency and the spatial 
domain. The properties of H and G in Mallat's pyramidal algorithms can be derived without 
reference to the multiresolution analysis [2]. 



Knowing the approximation function xm( 0 and all the detail functions di(t) at resolution 
level i = M, M+1, M+2, .... m-2, m-1, is necessary and sufficient to construct the 
approximation function Xm(t). The approximation of x(t) at resolution level m can be written as 



^m,n ^M,n 



neZ 



neZ 



M < Z < m neZ 



(3.8) 



for any M < m. Without loss of generality (WLOG), one may set M = 0 and ^(t) = (J»o,o(0 (the 
lowest resolution function physically realizable), and consider resolution level 0 < m < rnmax- In 
this report, we leave M as an arbitrary integer. An example pair for (j> and is the indicator 
function on [0, T) and the Haar wavelet on [0, T): 



V(0 = 



1/VT , 0 < t < T/2 
-1/VT , T/2 < t < T 



f l/>Tf 

m=\ 

0 



, 0 < t < T 
, otherwise 



(3.9) 



oo 



oo 
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In general, J <{)(t) dt = 4>(0) = VT, <1>(±«>) = 0, and J \|/(t) dt = 'F(O) = 0; the scaling 

-CX) -oo 

function <{)(t) can be interpreted as a lowpass function while \jr(t) is bandpass. The Haar wavelet 
forms the only orthonormal basis of compactly supported wavelets for which the associated 
scaling function ({) has a symmetry axis. Also note that the Haar wavelet is not continuous so a 
small error in the coefficients may cause a large fluctuation in the representation. Continuous and 
n times differentiable wavelets are known [2]. 

The two equivalences in (3.8) are valid for finite-energy deterministic functions in 
L2(IR); they also make sense for finite-power deterministic function in compact ICR, 

and for finite-energy and finite-power random processes w.p.l on compact ICR when compact 
scaling function <})(t) and compact wavelet function \jt(t) are employed [7], [12]. 

3.3 Wavelet Conditions and Generations 

It is clear that one needs to obtain the wavelet filter sequences h and g in order to perform 
the DWT. For the computation of the CWT in (2. 1 ), however, one needs to know the wavelet 
function \jr(t) explicitly. For compactly supported wavelets, knowing h and g can produce the 
time-functions <l)(t) and \|/(t). In this section, we discuss the conditions for generating an 
orthonormal scaling/wavelet pair, and give details for how to produce the time functions. This 
section is important for the fact that it provides working knowledge for both the DWT and the 
CWT. 

3.3.1 Wavelet conditions 

Consider the nested sequence of closed subspaces { Vm; m e 2 } as discussed in section 
3.1. The scaling function (at level m) (})m(t) is in Vm, which in turn is a proper subset of the 
closed space Vm+i. The basis of Vm+i is given by {<t>m+l,n(t); n e 2 } = t _ n); 

n e Z }, where we have set T = 1 WLOG. Since (t>m(0 is also in Vm+i, we can represent it in 
terms of the basis functions of Vm+i: 



4>m(t) — 4*m+l,n(0 > <}>m+l,n(0 

n 



(3.10) 



where the inner product < <{)m(t). 4>m+l,n(t) > is given by 



4*m(0> <l>m+l,n(0 ^ ~ J 4*m(0 <I>*m+l,n(0 dt 
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- j 2m/2<j)(2mt ) 2(m+l)/2,j,*(2m+lt _n) dt 

= V2 J<Kt)<K2t-n)dt (3.11) 

for real basic scaling function Equation (3. 1 1) describes exactly how two consecutive 

subspaces Vm C Vm+i are correlated. Note that (3.1 1) depends only on the scaling function <l)(t), 
but not on the individual levels m and m+1; it holds for any value of m. By defining h(n) via the 
following dilation equation: 

h(n) = J(})(t) (H2t -n) dt (3.12) 

for all n, we can rewrite (3.10) as 

<>m(t) = V2 Xh(n) 4>m+l,n(t) (3.13) 

n 

The dilation equation of (3.12) is the most fundamental equation in wavelet theory. One can 
proceed along this line of argument and derive the DWT of Figure 3.3. For a scaling function (|)(t) 
that has a finite support, the number of solutions to (3. 12) is finite, and so h is of finite length. 

It can be shown that [16]: 

|<Kt) dt = constant (3.14) 

and WLOG we can take the constant to be 1. This implies that <I>(f) = 1 at f = 0, which means that 
the scaling function is an impulse response of a lowpass type filter. Integrating both sides of 
(3.13) with respect to (w.r.t.) time t and letting m = 0 yields: 



1 = <2 Xh(n) J2l/2<t)(2t-n) dt 

n 



(3.15) 
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= 2 h(n) 



n 



Equation (3.15) is the first wavelet condition, which is a consequence of the nested nature of the 
increasing subspace sequence { Vn,}. For orthonormal scaling function <})(t), {<{)(t-nT); n e Z } is 
an orthonormal basis for Vq, and therefore, 



6(n) = <(|)(t),<|)(t-n)> 

= <2 2h(k)<)(2t-k), 2 X hG-2n) <K2t-j) > 
k j 

= 4$^ S h(k) h(j-2n) i 5(k-j) 
k j ^ 

= 2^ h(k)h(k-2n) 
k 



(3.16) 



Equation (3.16) is the second wavelet condition; it must hold if the scaling/wavelet system is 
orthonormal. Performing a discrete Fourier transform on both sides of (3.16) yields the conjugate 
filter equation of (3.5). Note that as a corollary to (3.16), we get by setting n = 0 that 
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(3.17) 



One needs only to find h such that it satisfies the two wavelet conditions of (3.15) and (3.16). 
Once h is known, the sequence g is given by g(n) = (-1)'“" h(l-n). These two filter sequences 
are all that are required to perform the DWT. 

3.3.2 Generation of wavelet function 



In the previous section, we have shown that if the scaling function <{)(t) is known, then we 
can determine h and g and perform the DWT with ease. Conversely, if h is found by solving the 
two wavelet conditions, then we can generate <{)(t) and \|i(t) recursively. The construction is as 
follows. Using (3.13) with m = 0, we get the fundamental recursive equation for the scaling 
function <})(t); 
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h(n) (|)(2t-n) 



(3.18) 
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Assuming that h = [h(0), h(l), h(2), h(3), h(M-l)], i.e., the h sequence has finite length, 

so that (})(t) has a compact support on [0, M-1], i.e., <})(t) is non-zero on this compact interval. We 
first find the values of <))(t) on the integers; this is done by using (3.18) to set up the following 
matrix equation; 
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(3.19b) 



for odd M. In general, solving (3.19a) or (3.19b) is equivalent to finding the eigenvector 
corresponding to the eigenvalue X = 1 of the matrix A given in below: 



^ = A ^ 



(3.20) 



We may impose the additional condition on ^ that <t)(n) = 1 . When the matrix A has zero 

n 

determinant there can be multiple X = 1 solutions to (3.20). Once we know the values of <})(t) on 
the integers, we can compute its values on the ^-integers by using the fundamental recursion of 

(3.18) as follows: 

4>(|) = 2 X h(n)<})(k-n) (3.21) 

n 

and so-forth for on the -j^, ^, . . . -integers, recursively. Thus, the scaling function <})(t), 

which has support on the compact time interval [0, M-1], can be determined explicitly for t being a 
dyadic number. Since the dyadic numbers (i.e., multiples of 2“*^, for all k > 0) are dense in E, the 
function <))(t) can be determined at all times t in the limit. 



2i: 

The wavelet function \|i(t) is similarly computed as follows. Since V|/(t) e Wq C Vj, we 
can represent t^/(t) as a linear combination of the basis functions of V i : j 

V(t) = Z < V(0. 21/2 4>(2t-n) > 21/2 <j)(2t-n) 
n 

= 2 X g(n)4>(2t-n) (3.22) 

n 

which is a recursion similar to (3.18). Since we already know (J»(t) on the dyadic numbers, it 
follows that \|T(t) can be computed using (3.22). The filter sequence g is a flip-altemate-negation 

of h, and we can set g = [g(0), g(l), g(2), , g(M-l)] so that V|f(t) also has a compact 

support on the time interval [0, M-1]. 1 

At this point we have shown how to generate \|/(t) for the CWT. The question is; How do 
we find h(n) in the first place? Although one may guess or even trial and error to find a h j 

sequence that satisfy the two wavelet conditions, in reality it is a difficult task to find a 
scaling/wavelet pair that has nice analytical properties. Nonetheless, several researchers have j 

succeeded in systematically generating wavelet functions. We show some of their work in the 
following. 

3.3.3 Compactly supported wavelets examples 



First of all, the simplest compactly supported wavelet is the Haar wavelet. It has a compact 
support on [0, 1], with h = [1/2, 1/2] and g = [1/2, -1/2]. The Haar scaling/wavelet pair is shown 
in Figure 3.4. a. Note that the functions are not continuous. 



Daubechies [2] has proposed a class of compactly supported orthonormal wavelets and she 
gave the hs for M = 4,6,8,10,12,14,16,18,20. For example, the Daubechies 4-coefficient wavelet 
has h = [h(0), h(l), h(2), h(3)] given by: 



h(0) 

h(2) 



1+V3 

8 

3-V3 

8 



h(l)= ^ 
h(3) = ^ 



(3.23) 



Note that h(0)+h(l)+h(2)+h(3) = 1, and h^(0)+h2(l)+h2(2)+h2(3) = 1/2. Substituting h into 
(3.19) we find the eigenvector for X = 1 to be 

<1)(1). <1)(2), <j)(3)] = [0, 1.366, -0.366, 0] (3.24) 

The corresponding scaling function <])(t) and the wavelet function \|/(t) (D4) are recursively 
computed, using (3.18) and (3.22), respectively, and they are plotted in Figure 3.4.b We also 
show the Daubechies 10-coefficient scaling/wavelet (Dio) pair in Figure 3.4.c and the 20- 
coefficient scaling/wavelet (D20) pair in Figure 3.4.d Observe that D20 looks like a modulated 
Gaussian. Also note that the Daubechies wavelets are increasingly smoother in time and better- 
contained in the frequency sp^ectrum. 

Pollen [19] has shown that it is possible to parametrize compactly supported wavelets. In 
particular, for M = 6, using a pinched torus, he was able to give a parametriztion of the h sequence 
as follows; 

h(0) = [(1 + cos a + sin a)(l - cos p - sin p) + 2 sin p cos a ]/8 
h(l) = [(1 - cos a + sin a)(l + cos P - sin P) - 2 sin P cos a ]/8 
h(2) = [1 + cos (a-p) + sin (a-p)]/4 
h(3) = [1 + cos (a-p) - sin (a-P)]/4 
h(4) = [l-h(0)-h(2)]/2 

h(5) = [1 -h(l)-h(3)]/2 (3.25) 

where < a, p < Tt are the two parametrizing angles. It is indeed a very nice and important 

result. As an example, the Haar wavelet has parametrization a = P, for any a. As another 

example, we randomly generate a = 0.2172 and P = 0.3359 and use (3.25) to obtain 

h = [0.0055, -0.0322, 0.4686, 0.5278, 0.0259, 0.0044]. (3.26) 

It turns out that this is a valid wavelet sequence; we carry out the recursive generation and we plot 
the resulting scaling function and wavelet function in Figure 3.4.e. This wavelet has dual 
passbands. Note that the scaling function and the wavelet function are orthonormal basis 
functions. One can experiment with (3.26) to generate very crazy looking wavelet functions. 
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Figure 3.4. a The Haar scaling function and wavelet function. 
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Figure 3.4.b Daubechies 4-coefficient (D 4 ) scaling function and wavelet function, 
and their energy spectra. 
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Figure 3.4.C Daubechies 10-coefficient (Dio) scaling function and wavelet function, 
and their energy spectra. 
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Figure 3.4.d Daubechies 20-coefficient (D 20 ) scaling function and wavelet function, 
and their energy spectra. 
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Figure 3.4.e A 6-coefficient scaling function and wavelet function with Pollen 
parametrization a = 0.2172, p = 0.3359, and their energy spectra. 
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3.4 Multiple-Phase Discrete Wavelet Transform 

In digital signal analysis, what one has at hand is a sequence of N samples, x(nT), where T 
is the sampling period and n = 0, 1 , 2, . . . Ns-1. It is convenient to consider an auxiliary function 
xo(t) in Vq: xo(t) = ^ bo,n 4>0.nO) by setting bo.n = x(nT) at resolution level 0, which is the 

n 

highest level available. Using a recursive filtering algorithm, the so-called discrete wavelet 
transform (DWT) of section 3.2, one can compute the lower-level resolution coefficients as (Fig. 
3.3): 



bm-l(i)= X>^2h(k) bm(2i + k) 
k 

dm-l(i) = X 2 g(k) bm(2i + k) (3.27) 

k 

for m = 0, -1, -2, -3 Here, h = [ h(-M+2), h(-M+3), . . . , h(0), h(l) ] and g = [ g(-M+2), 

g(-M+3), . . . , g(0), g(l) ] are the appropriate conjugate-quadrature discrete-time filter responses. 

oo 

The h coefficients are obtained by the dilation equation h(k) = 1 4>(t) 4>(2t - k) dt. We have set 

-OO 

g(k) = (-l)-M+3-k h(-M+3-k) so that the time indexes for the b’s and d’s are always non-negative. 
We also assumed that both ^ and \j/ have compact supports so that h and g have only finite number 
of filter taps ( = M). Roughly speaking, the approximation coefficients bm,n capture the 
approximate information in the signal sequence at resolution level m, while the difference 
coefficients dm,n retain the detail information going from resolution level m - 1 to m. The DWT 
has important applications in signal analysis, particularly in compression, detection, and 
classification. The choice of 4> and \j/ (hence h and g) is such that the signal is localized as much as 
possible in both time and frequency. Unfortunately, the DWT is not time invariant. If the signal 

sequence is time-shifted by one unit to become b = [ 0, bo.o. bo,i, bo,2. bo,N-l 1 then the 

corresponding DWT is not time-shifted and it will be completely different from the old DWT. This 
will not present a problem to the reconstruction of the signal sequence; however, the merits of 
DWT for signal detection and classification may be diminished. 

We note that for dyadic DWT, there are two different sets of DWT at each level; we name 
them as “phase-0” and “phase- 1” decompositions. Referring to Fig. 3.5, for phase-0 
decomposition, the 0-th coefficient at resolution level m - 1 is [h(l)] • [bm.ol. Ore 1-st coefficient is 
[h(-l), h(0), h(l)] • [bm,o, bm.i, bm,2) where denotes dot-product, the 2-nd coefficient is 



[h(2), h(-l), h(0), h(l)] • [bm.i, bm.2. bm,3. bm, 4 ], and so forth. On the other hand, the phase-1 
coefficients will read as [h(0), h(l)] • [bm.o, bm.i], [h(-2), h(-l), h(0), h(l)] • [bm.o, bm.i, bm.2. 
bm.sl. [b("2), h(-l), h(0), h(l)] • [bm,2. bm,3. bm,4. bm,5] . • • • • 
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Figure 3.5 Two phases of wavelet decomposition. 



We have proposed a multiple-phase DWT (MP/DWT) for signal analysis [10], [14], [13]. The key 
to this MP/DWT is to construct new filter sequences hm and gm (m < -1) at each decomposition 
level by inserting - 1 number of zeros between every two subsequent filter coefficients. 

The MP approximation coefficients at level m - 1, time indexed i, is given by the dot product of hm 



I I I 

and where ^ is the MP approximation coefficient sequence at level m denoted by ]^ = 

f t f 

K,i’ ^m.2 Specifically, we have 

1 

b ', . = 'V V 2 hm(k) b* . . , 

m-1,1 X j niv / m,i+k-l 



k= -M+Z-tM-lKZ-”-!) 



and 



m-l,i 



1 

^ V 2 gm(k) b^ j 

k= -M+2-(M-l)(2""-l) 



(3.28) 



The manner that the phase-0 and phase- 1 coefficients are interleaved at each level may be described 
by the tree diagram shown in Fig. 3.6. For a k-level decomposition, there are altogether 2^ 
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possible phases of decomposition. Each path in the tree of Fig. 3.6 represents an unique phase- 
sequence. One may choose to decompose the signal along any one of these 2 ^ phase-sequences 
and obtain a different DWT. Interleaving all of them yields our MP/DWT. Note that knowing the 
phase-sequence, one can trace back in the tree and reconstruct the signal from the DWT. Note that 
MP/DWT is a redundant transformation. It may be regarded as a discrete approximation of the 
continuous wavelet transform for xo(t). Its time-invariant property may be useful for some signal 
detection problems. The MP/DWT can be used to detect PSK signals, but we shall report the 
results elsewhere. For data compression applications, one can maximize the data compression 
ratio by optimizing over the 2*^ phase-sequences. 
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0 1 1 





Figure 3.6 Different phases of wavelet decomposition. 
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4 . DETECTION OF PSK SIGNALS 

We now use the CWT of (2.1) to analyze PSK signals. Viewing the CWT in a time-scale 
plane, we attempt to identify qualitative and quantitative features that would reveal the exact 
occasions of the phase shifts. First, we consider a basebad rectangular wave and its noisy 
versions. Next, we analyze a 180° phase shift in a sinusoid and its noisy versions. Haar, 
Daubechies, and complex modulated Gaussian wavelets will be employed in the decompositions. 

4.1 Rectangular Wave 

4.1.1 Haar wavelet 

A baseband rectangular wave and its CWTs using the Haar wavelet are shown in Figure 

4. 1 .a. The rectangular wave signal has a normalized duration of 1 second. The signal sequence 
x(nT) has length equal to 65 and it changes sign at h = 9, 25, 41, and 57, which corresponds to t = 
nT = 9/64, 25/64, 41/64 and 57/64 seconds, respectively. The Haar wavelet vjt(t) = \j^o(0 was 
shown in Figure 4.1.a. The transforms CWx(nT, 2"'") is computed using (2.1) for n = 0,1, . . . . 
. ,64, and m = 0, 1, 2, . ... ,6; they are shown in sequential order in Figure 4.1. a. Note that 
there are some end effects due to the finite signal window and we would ignore these end effects. 
At low level of decomposition, i.e., m = 0,1, or 2, the CWT has good frequency resolution and 
we see that the oscillating frequency of the signal is revealed. At higher level of decomposition, 
the time resolution is getting better and the CWT reveals the sign change instances with increasing 
accuracy. Specifically, at m = 6, a positive "spike" in the CWT shows a +1 to -1 sign change, 
and a negative "spike" shows a -1 to +1 sign change. The Haar wavelet is matched to the sign 
change so the good time resolution is expected. However, at high time resolution, the wavelet 
filter has large bandwidth. When there is noise, a lot of noise energy will be collected and the 
signal CWT will be smeared, as we will show in the next figure. 

We add noise to the rectangular wave and show the CWTs in Figure 4.1.b and 4.1.c. 
Power of the signal sequence is E(x^(nT)) = 1. The random noise sequence is Gaussian with 
mean zero and variance s^. Define the signal-power to noise-power ratio as: 

_ signal sequence power 
~ noise sequence power 



(4.1) 



In Figure 4.1.b, the SNR is 10 dB and in Figure 4.1.C the SNR is 5 dB. We see that the low level 
CWTs are not severely affected so we still have good frequency resolution. At higher level, the 
CWTs are increasingly smeared by the noise, especially at m = 6. This is because the output 
signal-power to noise-power ratio, SNRq, defined by: 



SNRo 



filtered signal sequence power 
filtered noise sequence power 



(4.2) 



is decreasing exponentially with m. Note that the wavelet filter (with impulse response V|/m(0) has 
bandwidth B = 2*2"' Oq and so the filtered noise power is increasing exponentially with m for 
wideband noise. The filtered signal power, however, will decrease with m since the r.m.s. center 
frequency of the wavelet filter is moving towards f = -h» and smaller and smaller amount of signal 
energy will be collected. If the signal is a wideband spread-spectrum signal instead of a 
narrowband signal, then it is possible that SNRq will remain large for many levels of 
decompositions. This suggests that WT is suitable for detecting spread-spectrum signals. 



For SNR = 1 0 dB as in Figure 4. 1 .b, one can observe from the CWT at m = 4 and m =6 
that there are 4 sign changes. However, when one attempts to pinpoint the exact times of 
occurrence by viewing the CWT at m = 6, only the changes at t = 9/64 and 41/64 can be located 
with confidence. 



For SNR = 5 dB in Figure 4. 1 .c, the noise power is too much for the CWT to reveal good 
time-resolutions. Only the oscillating frequency of the rectangular wave can be revealed. The 
Haar wavelet does not work well in noise because its poor specu-al properties, i.e., its spectrum 
decays too slow and has too many sidelobes. 
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Figure 4,1 .a CWT of rectangular wave with Haar wavelet. No noise. 
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Figure 4.1.b CWT of rectangular wave with Haar wavelet. SNR = 10 dB. 
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Figure 4. 1 .c CWT of rectangular wave with Haar wavelet. SNR = 5 dB. 



4.1.2 Daubechies wavelet 



Daubechies wavelets have better frequency properties than the Haar wavelet. However, 
they are not matched (in time) to the sign change in the rectangular wave. We show the CWTs for 
the rectangular wave with Daubechies 10-coefficient wavelet (Dio) in Figures 4.2. a, 4.2.b, and 
4.2.C for the noiseless, SNR = 10 dB, and SNR = 5 dB cases, respectively. The wavelet is first 
normalized in time so that it is confined to the interval [0, 1] at level m = 0. 

When there is no noise, the Daubechies wavelet is also capable in locating the sign changes 
with high time-resolution as shown in Figure 4.2.a. Furthermore, because its frequency spectrum 
is more confined than the Haar wavelet, the CWT at low level reveals the oscillating frequency of 
the rectangular wave with greater resolution. 

When SNR = 10 dB, the performance is worse, as shown in Figure 4.2.b. At levels m = 3 
and m = 4, it is possible to identify the first and the third sign changes correctly, just as in the Haar 
case. But there is a false indication around t = 0.4. When SNR = 5dB, this Daubechies wavelet 
is seen to perform not much better than the Haar wavelet, as shown in Figure 4.2.c. 

The problem with both Haar and Daubechies, as well as other unmodulated wavelets is that 
they are not tumable. It is difficult to select or "customize" the passband center frequency for a 
specific signal, thereby the wavelet filter may miss most of the signal energy. It is also not 
possible to have a relatively narrow bandwidth at high frequencies. The complex modulated 
Gaussian wavelet is capable to alleviate these shortcomings. 
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Figure 4.2.a CWT of rectangular wave with Daubechies 10-coefficient wavelet. 
No noise. 
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Figure 4.2.b CWT of rectangular wave with Daubechies lO-coefficient wavelet. 
SNR= 10 dB. 
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Figure 4.2.C CWT of rectangular wave with Daubechies 10-coefficient wavelet. 
SNR = 5 dB. 
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4.1.3 Complex modulated Gaussian wavelet 

By changing the modulating frequency, fv|/, of the complex modulated Gaussian wavelet, 
one can place the wavelet passband center frequency practically anywhere on the positive 
frequency axis. Changing fy, as well as a, in (2.10) gives us two extra dimensions to expose the 
signal features. Although there may be a lot of redundancy, this approach essentially gives us a 
tunable, variable-bandwidth multiresolutional analysis of the signal. It is also possible to adapt the 
analysis to different types of signal and noise. 

For simpler presentation, we fix a = 0.12 and only vary fy to be 4, 8, and 16 Hz. For the 
noiseless case, we show the corresponding time-scale decompositions (which are complex valued) 
in Figures 4. 3. a - 4.3. f. The first three figures depict the magnitudes of the CWTs, while the last 
three figures show the phases of the CWTs. Consider the magnitude plots in Figures 4.3.a, 4.3.b, 
and 4.3.C. There are altogether 3x7 = 21 snapshots of the wavelet filtered rectangular wave. Each 
set of 7 snapshots start with a different center frequency fy and the filter banks scan the frequency 
axis as m increases from 0 to 6. Most of these 21 magnitude plots do correctly reveal the sign 
change times (with different time-resolutions) and would be useful in the noisy cases. The phase 
plots in Figures 4.3. d, 4.3.e, and 4.3.f confirm the sign change times, and in addition they also 
reveal the polarity (+ to - or - to +) of each change. Note that when there is noise and the 
magnitude plots are not conclusive, we can consult the phase plots for a confirmation. This 
alternative source of information is not available in decompositions using real wavelets. 

For the case when SNR = 10 dB, we show the magnitude plots for fy = 4, 8, 16 in 
Figures 4.4. a, 4.4.b, and 4.4.c. In Figure 4.4.b at m = 6, the magnitude plot of the CWT clearly 
indicate that there are sign changes at the four correct instance, although it addition there is a 
possible false indication at t « 0.4. Figure 4.4. a at m = 3, 4, and 5 confirm that the findings are 
due to the signal but not the noise. Note that in Figure 4.4.C, the magnitude plots are not 
informative and certainly are not conclusive; the choice of fy = 16 is too large for this particular 
signal. 



When SNR = 5 dB, we show one set of magnitude plots of the CWT in Figure 4.5 with fy 
= 8. At level m = 6, the CWT reveals correctly the first, second, and the third sign change, 
probably misses the forth one and falsely add one at t » 0.4. The results can then be referenced 
with other magnitude plots using a different modulating frequency. Comparing with the Haar and 
the Daubechies wavelets, the performance of the complex modulated Gaussian is better, especially 
in the noisy cases. 
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Figure 4.3.a Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, f\j/ = 4. No noise. 
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Figure 4.3.b Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet. o = 0.12, fy = 8. No noise. 
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Figure 4.3.C Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, f^j, = 16. No noise. 
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Figure 4.3. d Phase of CWT of rectangular wave with complex modulated 
Gaussian wavelet, a = 0.12, fvj, = 4. No noise. 
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Figure 4,3.e Phase of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, f,j/ = 8. No noise. 
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Figure 4.3.f Phase of CWT of rectangular wave with complex modulated 
Gaussian wavelet, a = 0.12, f\j/ = 16. No noise. 
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Figure 4.4.a Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, fy = 4. SNR = 10 dB. 
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Figure 4.4.b Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, f\)/ = 8. SNR=10dB. 
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Figure 4.4.C Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, a = 0.12, fvjf = 16. SNR = 10 dB. 
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Figure 4.5 Magnitude of CWT of rectangular wave with complex modulated 
Gaussian wavelet, o = 0.12, fv}/ = 8. SNR = 5dB. 
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4.2 180* phase shift 

The signal to be analyzed is now changed to a sinusoid with a 180' phase shift at t = 0, as 
shown in the fu-st subplot of Figure 4.6. We use the Haar, Daubechies lO-coefficient, and the 
complex modulated Gaussian wavelets to analyze this signal and plot their CWTs in Figures 4.6, 
4.7, and 4.8, respectively. For the Haar wavelet, the phase change also app>ears as a phase change 
in most of the CWTs as shown in Figure 4.6. The Daubechies wavelet does a better job in 
identifying the phase change, as shown in Figure 4.7. The complex modulated Gaussian wavelet 
also reveals clearly the phase shift in the magnitude plots of the CWT at m = 4 and m =5, in 
addition a phase change is confirmed by m = 1 and m =2 plots, as shown in Figure 4.8. 

When there is noise., we found that both the Haar and Daubechies wavelets failed to 
perform. 

When SNR = 10 dB (note that the signal power is equal to 0.5), we choose a = 0.12 and 
f\|, = 12 for the complex modulated Gaussian and show the magnitude and phase plots in Figures 

4.9. a and 4.9.b. In Figure 4.9. a, the magnitude plot of the CWT at m = 2 reveals a phase shift 
activity somewhere on or before t = 0.5. (Note that there is an end-effect in the DWT due to finite 
signal window.) But the rest of the magnitude plots fail to pinpoint a better time location. When 
we consult the phase plots as shown in Figure 4.9.b, the m = 6 plot suggests that there is an 
abrupt change of phase at t = 0.5. This is confirmed by the m = 3 and m = 6 phase plots in Figure 

4.9. c, where we use a different fy = 16. 

When SNR = 5 dB, we choose a = 0. 12, fy = 8. The magnitude subplot for m = 2 of 
Figure 4.10.a reveals some phase activities; however, the other magnitude subplots and the phase 
plots of Figure 4.10.b yield no confirmative information. Detection of a single phase shift is very 
sensitive to the noise. However, when the signal window is much longer and more phase shifts 
are present, the complex modulated Gaussian would perform due to a larger output signal-power to 
noise-power ratio (SNRo, see (4.2) ). 
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Figure 4.6 Magnitude of CWT of 180' phase shifted sinusoid with Haar 
wavelet. No noise. 
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Figure 4.7 Magnitude of CWT of 1 80’ phase shifted sinusoid with Daubechies 
lO^oefficient wavelet. No noise. 
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Figure 4.8 Magnitude of CWT of 1 80* phase shifted sinusoid with complex 
modulated Gaussian wavelet, a = 0. 12, fy = 8. No noise. 
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Figure 4.9.a Magnitude of CWT of 1 80' phase shifted sinusoid with complex 
modulated Gaussian wavelet, o = 0.12, f\(( = 12. SNR = 10 dB. 
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Figure 4.9.b Phase of CWT of 1 80“ phase shifted sinusoid with complex 

modulated Gaussian wavelet, a = 0.12, fy = 12. SNR = 10 dB. 
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Figure 4.9.C Phase of CWT of 180' phase shifted sinusoid with complex 

modulated Gaussian wavelet, o = 0.12, fy = 16. SNR = 10 dB. 
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Figure 4.10.a Magnitude of CWT of 180' phase shifted sinusoid with complex 
modulated Gaussian wavelet, o = 0.12, f,j, = 8. SNR = 5 dB. 
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Figure 4. 1 0.b Phase of CWT of 1 80' phase shifted sinusoid with complex 

modulated Gaussian wavelet, a = 0.12, f\jf = 8. SNR = 5dB. 
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5. CONCLUDING REMARKS 

Wavelet decompositions have better time-resolution at higher frequencies. They have 
better frequency-resolution at lower frequencies. Wavelet decompositions are suitable for 
analyzing signals that have high-frequency short pulses and/or low-frequency long pulses. The 
discrete wavelet transform (DWT) is primarily useful for data compaction, although we have found 
some of its potentials in signal detection. The continuous wavelet transform (CWT) can be used in 
signal detection. We have shown that by using complex modulated Gaussian wavelets, the time 
and frequency features of a signal can be revealed via a tunable, variable-bandwidth, 
multiresolution analysis. We have shown that CWT can be used in detecting PSK signals. 
However, the performance is sensitive to the level of noise. 
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